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Abstract. We present a variational study of the helicity moduli of an anisotropic 
quantum three-dimensional (3D) XY-model of YBCO in superconducting state. It is 
found that both the a6-plane and the c-axis helicity moduli, which are proportional 
to the inverse square of the corresponding magnetic field penetration depth, vary with 
temperature T as T 4 in the zero temperature limit. Moreover, the c-axis helicity 
modulus drops with temperature much faster than the a6-planc helicity modulus 
because of the weaker Josephson couplings along the c-axis compared to those along 
the a6-plane. These findings are in disagreement with the experiments on high quality 
samples of YBCO. 
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1. Introduction 

In a recent paper [1] , following the suggestions of Roddick and Stroud [2] and of Emery 
and Kivelson [3-6], we have examined the possibility that a classical XY-model might 
explain both the aft-plane and the c-axis electromagnetic response of optimaly doped 
single crystal YBCO [7]. A Monte Carlo study was applied to a bilayer model for 
YBCO in which all the microscopic degrees of freedom are assumed to be integrated 
out except for the phase of the superconducting order parameter and all of the phase 
dynamics results from Josephson couplings within the layers and between the layers 
grouped in a stack of bilayers. We calculated the temperature (T) dependence of the 
helicity moduli which in this model correspond to the inverse square of the magnetic 
field penetration depths X a f,(T), A C (T), i.e. to the superfluid densities n S;a b(T), n sc (T). 
Our results were in sharp contrast with experiment [7]: while we found that both ab- 
plane and the c-axis superfluid densities decrease linearly with T at low temperatures, 
experimentally only the aft-plane superfluid density drops linearly with temperature at 
low T. The observed c-axis penetration depth A C (T) never has the linear temperature 
dependence found in aft-plane. 

Here we extend the model of Ref. [1] to include the quantum effects associated 
with the Coulomb charging energy due to the excess number of Cooper pairs within 
coarse-grained regions with linear dimensions on the order of the superconducting 
coherence length. As we are primarily interested in the low temperature helicity 
moduli of this anisotropic, but three-dimensional, model we apply the variational self- 
consistent phonon approximation of Wood and Stroud [8]. The inhomogeneity of 
Josephson couplings and layer spacings perpendicular to the layers makes the calculation 
of the z-axis helicity modulus a nontrivial extension of the work by Roddick and 
Stroud [2]. Moreover we specifically address the analytic form (power law behavior) 
of the temperature dependence in the in-plane and out-of-plane superfluid densities by 
calculating the relevant weighted phase-phonon densities of states. 

The rest of the paper is organized as follows. In section 2 we discuss our model 
and derive the expressions for the helicity moduli within the self-consistent phonon 
approximation. Section 3 contains our numerical results for the helicity moduli, an 
analysis of their low temperature behavior including the relevant analytical results for 
weighted phase-phonon densities of states and our conclusions. 

2. The model and helicity moduli 

We consider an anisotropic quantum 3D XY-model for a system with bilayer structure 
described by the Hamiltonian 



H =T + V ab + V c 



(1) 




d 



(2) 



d(j)i,i,s 
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V<*= S - cos(0 M , s - (f)j,i, s ))} , (3) 

l,s,{i,j) 

Vc = J2 ~ cos (0i,«,2 - &,J,l)) + ^1(1 - 008(^^+1,1 - 0i,i,2))] • (4) 

i,i 

Here, the sum over / runs over a stack of bilayers, the sum over s — 1,2 runs over 
two layers in a given bilayer, denotes the nearest neighbors within a single layer, 
the sum over % runs over the sites in a given layer and is the phase of the order 
parameter on site % of the layer s in the bilayer /. T is the charging energy associated 
with excess number of Cooper pairs on site (i,l,s) [8], while V a b and V c describe the 
Josephson coupling energy within the layers and between the layers, respectively. The 
model considered here differs from the one examined by Roddick and Stroud [2] in 
that the Josephson coupling constant between two layers within a given bilayer J± and 
their spacing q, are different from the Josephson coupling constant between layers in two 
adjacent bilayers J' L and the bilayer spacing d . As a result, the calculation of the helicity 
modulus along direction perpendicular to the bilayers is a nontrivial generalization 
of the work by Roddick and Stroud [2]. 

In this work we apply the self-consistent phonon approximation [8] to the low 
temperature thermodynamics implied by the Hamiltonian H. One would expect that 
for a three-dimensional model this approach provides a reasonable description of the low 
temperature thermodynamics, at least for not too large values of U / max{ J 1; Jj_, J ; ± }. 
However, one cannot a priori rule out a possibility that quantum Monte Carlo treatment 
reveals, in particular in the limit of high anisotropy, a qualitatively different behavior of 
the type found by Jacobs et al. [9] in purely two-dimensional case. We leave quantum 
Monte Carlo treatment of our model for the future investigations. 

In variational self-consistent phonon approach the trial free energy has the form 

F t = F h -(H- H h ) h , (5) 



where 



H h =T + V£ ) + VW, (6) 

v£ fc) = E ^(^s-ha, (7) 

l,s,{i,j) 



7, Wi,l,2 - 9i,l,l) + 7, C0i,I+i,i ~ 0i,Z,2) 



(8) 



l,i 

(■ • -)h denotes the statistical average with respect to the quasi-harmonic Hamiltonian 
Hh with temperature dependent "spring constants" , and is the Helmholtz free energy 
defined by the Hamiltonian Hh- The values of Ki(T), K±(T) and K' ± (T) are determined 
from self-consistency equations which result from minimizing the trial free energy F t . 

Transforming to the normal modes and using (cos(0 — <f>'))h = -Re(exp(i(0 — 4>')))h 
together with Mermin's theorem [10] we obtain 



F t = k B T ln f 2sinh ^l+ 4iV 



s=l,2 k 



J 1 (l-e-^)-^K 1 D 1 
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+N 



) - - 2 K ± D ± 



+ 



Ji(l-e-^)_I^l]} 



and 



■Mr) 



K^T) = J ie -t J 

k' ± (t) = j; e -^i( T » , 



(9) 

(10) 

(11) 
(12) 



where 








Di(T) = 


E E (sin 2 (fc x a/2) + sin 2 (A; 2/ a/2 

iV s=l,2 k 


)) h coth^ s(k) 
') 2Mu s (k) 2ksT ' 


(13) 


D ± (T) = 


i \ - \ - h 







X 



1 + (-!)< 



X_l + K' ± cos(k z c) 



hu s (k) 

J(K ± + K' ± ) 2 - 4K ± K' ± sin 2 {k z c/2) ~ ' 2/c b T 



coth 



(14) 



N 



EE 



ft 



- , , k 2Mu 8 (k) 
l + (-l) 



x 



X_l cos(£; z c) + K' ± 



M 



(K ± + K' ± ) 2 - 4K ± K' ± sin 2 (k z c/2) 
K ± + K> ± 



: coth 



2k B T 



(15) 



2K X (sm 2 (k x a/2) + sva 2 (k y a/2)) + 



: yJ(K ± + K' ± y - 4K ± K' ± sin 2 (k z c/2) 



1/2 



(16) 



The sums over k run over the first Brillouin zone, N is the number of unit cells and 
M = h 2 /(AU). Equations (9-16) can be solved iterativly at fixed temperature for a 
given set of the bare coupling constants J±, J±, J' ± and for a given Coulomb parameter 
U. In (16) index 1 denotes the acoustic phase phonon dispersion and index 2 denotes 
the optic one. 

When a uniform vector potential A is applied its effect on the Hamiltonian is to 
shift the phase difference between points r x and r 2 by A 1>2 = 2nA ■ (r 2 — ri)/$ , where 
$o — hc/2e is the flux quantum. For example, if A is perpendicular to the bilayers 
&,j,2-0i,j,i and <fri,i+i,i-(j)i,i,2 in Eqs. (4) and (8) are replaced by <j>i, l , 2 -<i>ij,i + (2n/&o)c b A 
and <pi,i+i,i — 4>i,i,2 + (27r/$o)c / A, respectively. As a result, for A perpendicular to the 
bilayers the quasi-harmonic Hamiltonian, (6-8), takes the form 



H h (A) = H h (0) + H'(A) 
2tt 



0k=O,l, 



(17) 
(18) 



where a constant term, proportional to A 2 , has been dropped and 0k,s is the Fourier 
transform of the phase variable 0i iS defined on Bravais lattice sites 1 (s=l,2). 

The computation of the helicity modulus along direction perpendicular to the 
bilayers requires even more care than in the classical case [1] because dHh(A)/dA 
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does not commute with Hh{A) (see Eqs. (17-18)) and therefore <9exp(— (3Hh{A)) / dA ^ 
-(3(dH h (A)/dA)exp(-(3H h (A)). Hence the helicity modulus d(J(A)) hA /dA\ A=0 , 
where J (A) is the current operator through a particular bond perpendicular to 
the layers and (■■■)h,A is the usual quantum statistical average given by Hh(A), 
is not given by the standard fluctuation-type formula {dJ{A) / dA\A=o)h,A=o — 
{l/k B T)(J{Q)dH h {A)/dA\ A=Q ) hA=0 + {l/k B T)(J{^ 

Instead, we compute the average current {J(A)) h ^o to the lowest order in A using 
the perturbation theory and evaluate the helicity modulus from d(J(A))h,A^o/ dA\ A =o- 
Focusing on a bond between two layers in a bilayer (strong bond) we have [1] 

J S ( A ) = J± sm(0M, 2 - + (2vr/$ )c fe A). (19) 

Writing exp(— (3Hh{A)) = exp(— j3Hh(fyU(/3, A) we get the equation of motion 
dU((3,A)/d(3 = H'((3,A)U((3,A), H'(f3,A) = exp(f3H h (0))H'(A) exp(-/3H h (0)), and 
the boundary condition U(0, A) = 1. Thus, to the lowest order in A we have 

U(P, A) = 1 - / drH'(T, A)U(t, A) « 1 - / drH'(T, A) . 
Jo Jo 

Using this result, (sin(0 — <f>'))h,A^o = (l/i)Im(exp(i(<f> — 4>')))h,A^o and the operator 

identities 

£ e c(S+8t) = I9 c{ah+ht) + C eC(Mt) 

coa 2 
ftt e c(S+St) = 1 A e c(S+aSt)| _ £ c (S+8t) (21) 

cda 1 2 v ; 

valid for any two operators b and such that [b, W] — 1, we find to the lowest order in 
A 

'^' ^™ (l>+4 (22) 

Thus the z-axis helicity modulus of our model in quasi-harmonic approximation is given 
by 

lzziT) = %K ± (T) + K' ± (T) {Cb + C) - (23) 

Needless to say, we get the same result for 7 ZZ (T) by computing the average current 
along the weak bond, i.e. the one between the bilayers, as expected from Kirchhoff's 
first law (note the symmetry of the expression (23) under the exchange of (K±,Cb) and 
(K' ±: c% 

Another way to write our central result Eq. (23), which is intuitively appealing, 
would be 

lzz (T) = ^K ± , eff (T)(c b + c>), (24) 



where 



K±,eff K t K\ 



1 +^r- (25) 
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Equation (25) gives the effective spring constant for two linear springs with constants K± 
and K' ± which are connected in series. The lattice parameter in direction perpendicular 
to the layers is q, + d. The form given by Eqs. (24-25) is analogous to what one finds 
using quasi-harmonic approximation for a lattice without basis [2], e.g. for the in-plane 
helicity modulus of our model 

2n 

lxx (T) = —K l (T)a. (26) 
<±>o 

In that case the helicity modulus is essentially the stiffness with respect to the phase 
twist on neighboring sites. We note that within the quasi-harmonic approximation 
7 2Z (T) does not depend on the relative size of the bond lengths c b and d (see Eqs. (10- 
16), (23)), and the z- axis lattice parameter c = c& + d enters only implicitly as k z c in 
the Brillouin zone sums. 

In conclusion of this section we note that if one takes first the classical limit h — > 
and then the zero temperature limit T — > in Eqs. (10-16), (23) one recovers the classical 
results obtained previously for the bilayer model using the low temperature spin-wave 
expansion [1]. 



3. Numerical results and conclusions 



We have computed the helicity moduli 7 22 (T) and 7 XX (T) for different choices of the 
bare coupling constants Ji, J±, J' ± and since they give qualitatively very similar results 
we present here only the data for a single set of parameters J\ = 1, J± = 0.1, and J' ± 
= 0.001. The value of Josephson coupling between the layers of a bilayer is taken to 
be one order of magnitude smaller than the in-plane Josephson coupling and the value 
of Josephson coupling between bilayers J' ± is two orders of magnitude smaller than J±. 
We have considered five values of the Coulomb repulsion parameter U = 10~ 6 , 10~ 3 , 
10~ 2 , 10 _1 and 1 (in units of the a6-plane coupling Ji). Since the helicity modulus 
is proportional to inverse square of the magnetic field penetration depth A(T) [2], we 
present our results for the helicity moduli in Fig. 1 as A 2 (0)/A 2 (T) as a function of T. 
We did not attempt to determine the transition temperature T c for various values of U 
as it is unlikely that the self-consistent phonon approximation treats correctly vortex- 
antivortex pairs close to T c and we are primarily interested in the low temperature 
dependence of the helicity moduli. 

We point out that in the zero temperature limit 7 2Z (T) and ^ X x(T) are not linear 
in T for any finite U. It is easy to show that in the limit T — > both helicity moduli 
vary as T 4 . Indeed, the expressions for -Di(T), Dj_(T) and D' ± (T), which determine the 
temperature dependences of 7 ZX (T) and 7 2Z (T), can be conveniently written in terms of 
various weighted phase-phonon densities of states 

ZMT) = fi^M + FfM^ccth^ , (27) 
ff V) = T^\sJ FBZ d3k (sin 2 (M/2) + sin 2 (*:„a/2)) S m (u, - fiu,(k)) , s = 1,2 , (28) 
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T/J 1 



Figure 1. The temperature dependence of the helicity moduli for different values of 
the Coulomb repulsion parameter U. The top group of five curves are the a6-plane 
helicity moduli for different values of U, and the bottom group of curves are the 
corresponding c-axis helicity moduli. 
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D±(T) = f~MF?\u>) + F?\u)] — coth^ f , (29) 



o 



(2n) 6 Jfbz 



y/(K ± + K' ± Y - 4K ± K> ± sm 2 (k z c/2)_ 
x5 (3) (uj -fno s (k)) , 8 = 1,2, (30) 

D'AT) = r^ d ^ F ± ) ^ + F T\^)}^^th^, (31) 



o 



± V ; (2ti) 3 Jfbz 



x y K ± cos(k z c) + K' ± 



J{K ± + K' ± y - 4K ± K' L sin 2 (k z c/2)_ 
x5 (3) (uj -huj s (k)) , s = 1,2, (32) 

where v = a 2 c is the volume of the unit cell. 

We have computed Fi s \ F[ s) , F'} s \ s =1,2, numerically using the tetrahedron 
method [11]. At low temperatures it is the frequency dependence of the acoustic 
phase phonon weighted densities of states F^\ F^ and F^ that determines the 
temperature dependence of D\(T), D±(T), D' ± (T), and thereby of K\(T), K±(T), 
K'j_(T) (Eqs. (10-12)), which in turn determine the temperature dependences of the 
helicity moduli (Eqs. (23) and (26)). In the limit uu — > we find that because of the 
weighing factors F^\ F^ and F^ are all proportional to uj 4 (the coefficients are 
ai = J2(K~ ± + K{)/(UK ± K' ± )/2 7 n 2 3 (UK,) 2 , a ± = a 1 (2(K 1 /(K ± + K' ± )(K' ± /K ± ))), 
and a' ± = a±(K±/ K' ± ) 2 , respectively). Then, it is easy to see, by rescaling the integration 
variables in Eqs. (27), (29), (31), that D^T) - 0,(0) oc T 4 , etc.. After expanding, one 
finds that both helicity moduli are decreasing with T as T 4 in the zero temperature limit. 
We note that Xiang and Wheatley [12] found Independence of the c-axis superfluid 
density, which is proportional to j zz (T), due to entirely different physical mechanism 
than the one considered in the present work. 

The results shown in Fig. 1 are qualitatively different from what is found 
experimentally on single crystals of bilayer compound of optimally doped YBCO [7] 
(see Fig. 2 in [7]). While experimentally the afe-plane helicity moduli are linear in 
T at low temperatures we obtain Independence in the limit T — > for any finite 
U. Roddick and Stroud [2] found (for a one-layer model of YBCO) that including 
dissipation through coupling to an Ohmic heat bath restores the linearity of helicity 
moduli (presumably in all directions). We did not consider such coupling in this work. 
Physically the dissipation results from quasiparticle degrees of freedom [13], and we have 
discovered recently [14] that the effect of nodal quasiparticles on phase fluctuations in a 
d-wave superconductor might render the XY-model type of description of such system 
meaningles in the limit T — > 0. Nevertheless, the biggest discrepancy between the 
experiments [7] and our results in Fig. 1 is that the calculated r ) Z z(T)/^ zz (0) decreases 
much faster with temperature than the calculated j xx (T) /j X x(ty- This feature was also 
present in our Monte Carlo study of the classical version (U = 0) of the same model 
[1] and in either case it is a direct consequence of the weaker Josephson couplings 
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in direction perpendicular to the layers compared to the in-plane Josephson coupling. 
Dissipation, treated as in [2], will not affect that result. 

In conclusion, neither a classical nor a quantum-mechanical XY-model (at least in 
the self-cosistent phonon approximation) can consistently account for both afr-plane and 
c-axis electrodynamics of YBCO observed in experiments [7]. 
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